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ABSTRACT 

We present a model aimed to reproduce the observed spectral energy distri- 
bution (SED) as well as the ammonia line emission of the G31.41+0.31 hot core. 
The hot core is modeled as an infalling envelope onto a massive star that is un- 
dergoing an intense accretion phase. We assume an envelope with a density and 
velocity structure resulting from the dynamical collapse of a singular logatropic 
sphere. The stellar and envelope physical properties are determined by fitting 
the observed SED. From these physical conditions, the emerging ammonia line 
emission is calculated and compared with subarcsecond resolution VLA data of 
the (4,4) transition taken from the literature. The only free parameter in this 
line fitting is the ammonia abundance. The observed intensities of the main and 
satellite ammonia (4,4) lines and their spatial distribution can be well reproduced 
provided it is taken into account the steep increase of the gas-phase ammonia 
abundance in the hotter (> 100 K), inner regions of the core produced by the 
sublimation of icy mantles where ammonia molecules are trapped. The model 
predictions for the (2,2), (4,4), and (5,5) transitions, obtained with the same set 
of parameters, are also in reasonably agreement, given the observational uncer- 
tainties, with the single-dish spectra of the region available in the literature. The 
best fit is obtained for a model with a central star of ~ 25 Mq, a mass accretion 
rate of ~ 3 x 10^^ Mq yr~^, and a total luminosity of ~ 2 X 10^ Lq. The outer 
radius of the envelope is 30,000 AU, where kinetic temperatures as high as ~ 40 
K are reached. The gas-phase ammonia abundance ranges from ~ 2 x 10^^ in 
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the outer region to ~ 3 x 10"^ in the inner region. To our knowledge, this is 
the first time that the dust and molecular line data of a hot molecular core, in- 
cluding subarcsecond resolution data that spatially resolve the structure of the 
core, have been simultaneously explained by a detailed, physically self-consistent 
model. This modeling shows that hot, massive protostars are able to excite high 
excitation ammonia transitions up to the outer edge (~ 30,000 AU) of the large 
scale infaUing envelopes. 

Subject headings: Circumstellar matter — ISM: individual (G31. 41-1-0. 31), molecules 
— radiative transfer — stars: formation 



1. Introduction 

Hot molecular cores (hereafter HMCs) are small, dense, hot, and dark molecular clumps, 
usually found in the proximity of ultracompact HIT regions (e.g., Kurtz et al. 2000). Unlike 
HIT regions, these objects present weak or undetectable free-free emission. The lack of free- 
free emission has been interpreted as due to an intense mass accretion phase that quenches 
the development of an HII region (e.g., Walmsley 1995). Therefore, these HMCs may be 
the precursors of Hll regions, tracing the earliest observable stage in the life of massive 
stars. A detailed modeling of the emission of the HMCs is a necessary requirement to assess 
their physical nature and their relationship with the birth of massive protostars. Osorio, 
Lizano, & D'Alessio (1999) showed that the spectral energy distribution (SED) of the dust 
continuum emission of several HMCs without detectable free-free emission is consistent with 
that of a massive envelope collapsing onto a B type star, with an age less than 10^ yr, and 
accreting mass at a rate of 10~^-10~^ Mq yr"-*^. These intense mass accretion rates are 
high enough to prevent the development of an ionized region around the central star. The 
physical properties of the HMCs were inferred using the collapse of the Singular Logatropic 
Sphere (SLS; McLaughlin & Pudritz 1997) as a dynamical model for the envelope. The SLS 
collapse models were chosen because these envelopes are massive enough to reproduce the 
strong millimeter emission observed, without requiring excessively high mass accretion rates, 
thus preventing excessive emission in the near and mid infrared (see more details in Osorio 
et al. 1999). 

One can further test Osorio's models, including the velocity field of the accreting en- 
velopes, by comparing their expected molecular emission with observations of appropriate 
molecular tracers. High-excitation ammonia transitions are well suited for this kind of com- 
parison, since they trace the hot and warm dense gas in the HMCs, and can be observed 
with subarcsecond angular resolution using the Very Large Array ( VLA) . Such observations 
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are among the highest angular resolution molecular hne data of HMCs that can be obtained 
at present, allowing to determine the variations of the line intensity across the sources. Fur- 
thermore, the ammonia lines have hyperfine structure that is sensitive to the optical depth 
and permit an accurate measurement of the column density. 

In order to do this test, we selected the HMC located ~ 5" to the southwest of the peak 
of the G31.41+0.31 ultracompact HII region. Hereafter, we will call this source G31 HMC. 
This HMC is located at a distance of 7.9 kpc (Cesaroni et al. 1998). It is an interesting 
HMC because it exhibits strong millimeter continuum emission, as well as molecular emission 
of high excitation lines; it is associated with a group of water mascrs, and with a bipolar 
outflow (Cesaroni et al. 1994, 1998; Gibb, Wyrowski & Mundy 2004). The presence of high 
excitation molecular emission in this luminous source indicates that G31 HMC (together 
with a few other objects; Beuther & Walsh 2008, Longmore et al. 2007) is one of the 
hottest molecular cores discovered so far. The lack of strong free-free centimeter emission 
(Araya et al. 2003) indicates that G31 HMC could be in a phase of intense mass accretion, 
being unable to develop a detectable HII region. Therefore, we consider G31 HMC as a 
good candidate to test the basic hypothesis of Osorio's model, namely, that young massive 
stars are formed inside HMCs with very high accretion rates. In this model the onset of a 
detectable ultracompact HII region is quenched by the very high mass accretion rate. 

G31 HMC has a database of continuum flux density measurements covering a wide range 
of wavelengths although, unfortunately, most of the data points were obtained with relatively 
poor angular resolution and the SED is poorly constrained for A < 850 fira. G31 HMC is 
associated with line transitions from a variety of molecular species, such as HCO''", SiO 
(Maxia et al. 2001), ^^^CO (Olmi et al. 1996), CS (Anglada et al. 1996), H2S, C^^O (Gibb, 
Mundy & Wyrowski 2004), and CH3CN (Beltran et al. 2005). Particularly interesting is the 
association of G31 HMC with strong ammonia emission (Churchwell et al. 1990, Cesaroni 
et al. 1992) whose (4,4) inversion transition has been observed with subarcsecond angular 
resolution using the VLA (Cesaroni et al. 1998; see their Figs. 2c, 4c, and 9c). This makes 
G31 HMC one of the few sources where a high signal-to- noise ratio analysis of the spatial 
variation of the ammonia emission along the core can be carried out. The high intensity of 
the observed ammonia (4,4) emission, the ratio of the main to satellite lines close to unity, as 
well as the unusually large line widths indicate extreme physical conditions that suggest that 
this core is harboring an O type star. Because of this, we chose G31 HMC as a candidate 
to test the molecular emission predictions of the Osorio's model for the early stages of the 
formation of massive stars. 

We note that VLA observations at 7 mm (Araya et al. 2003) reveal the presence of a 
binary system at the center of G31 HMC. However, the model of Osorio et al. (1999) is 
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spherically symmetric, with a single central source of luminosity. This assumption is valid 
as long as the binary separation (~ 0''2) is much smaller than the observed size of the HMC 
(~ 8"). Likewise, deviations from spherical symmetry either by rotation (disks) or by an 
intrinsic elongation of the cloud are not considered by Osorio's model. To test non-spherical 
models, one must include high angular resolution mid-IR data that would allow to properly 
constrain parameters such as the degree of flattening of the envelope, the centrifugal radius, 
or the effect of a disk and/or cavities caused by outflows. This kind of data has been obtained 
for other HMCs and is highly sensitive to the geometry of the cloud (e.g., De Buizer, Osorio 
& Calvet 2005). Unfortunately, such observations are not yet available for G31 HMC. 

Observations of several molecular tracers reveal, however, the existence of small veloc- 
ity gradients in G31 HMC that have been interpreted as suggestive of either outflow (Gibb, 
Wyrowski, & Mundy 2004, Araya et al. 2008) or rotation motions (Beltran et al. 2005). 
Sensitive high angular resolution observations are needed to clarify this issue. Despite the 
possible presence of small velocity gradients, we expect the kinematics of the core to be 
dominated by the overall infall motions of the envelope. In this work we assume that the 
spherically symmetric models of Osorio et al. (1999) are adequate to study the dust contin- 
uum and molecular line emission of G31 HMC. 

The paper is organized as follows. In §2, we describe the modeling procedure, in §3 
we present the observational data available from the literature, in §4 we model the dust 
and ammonia data and we compare our results with the observations, as well as with other 
modeling attempts of this source reported in the literature. Finally, in §5 we summarize our 
conclusions. 



2. Modeling 

2.1. General Procedure 

Following Osorio et al. (1999), we model a HMC as a spherically symmetric envelope 
of dust and gas that is collapsing onto a recently formed central massive star. An accretion 
shock is formed at the stellar surface. Thus, the stellar radiation and the accretion luminosity 
provide the source of heating of the envelope. The temperature of the envelope is a function 

of the distance to the center. The density, infall velocity, and turbulent velocity dispersion 
for each point of the envelope arc obtained from the solution of the dynamical collapse of 
the SLS (McLaughlin & Pudritz 1997). The parameters of the SLS model arc chosen by 
calculating the emerging dust emission, assuming a constant dust to gas mass ratio along 
the envelope, and fitting the observed SED. 
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Once obtained the physical conditions in the envelope from the fit to the observed SED, 
the excitation of the ammonia molecule is calculated and the radiative transfer is performed 
in order to obtain the emerging ammonia spectra for different lines of sight towards the enve- 
lope. In these calculations, the gas-phase ammonia abundance, that is a function of radius, 
is the only free parameter. Wc consider two simple cases: uniform gas-phase abundance 
and uniform total ammonia abundance. In the latter case, the ratio of gas to solid phase 
ammonia as a function of radius is determined by calculating the balance between conden- 
sation and sublimation of molecules on dust grains. The gas-phase ammonia abundance is 
constrained by comparing the ammonia model spectra with the observed spectra towards 
different positions. 



2.2. Physical Structure of the Envelope 

For the collapsing envelope, we adopt the density, psLs(^)) infall velocity, fsLs(^); and 
turbulent velocity dispersion, a"sLs('")) distributions resulting from the solution of the self- 
similar collapse of the SLS (Adams, Lizano & Shu 1987, unpubhshed notes) . The SLS has a 
logatropic pressure, P — Poln(psLs/po), where Pq is a constant that sets the pressure scale, 
and po is an arbitrary reference density, introduced by Lizano & Shu (1989) to empirically 
take into account the observed turbulent motions in molecular clouds. The collapse of the 
SLS has been studied by McLaughlin & Pudritz (1996, 1997; hereafter MP97), and it has 
been generalized by Reid, Pudritz & Wadsley (2002), and Sigalotti, de Felice & Sira (2002) 
to include three-dimensional hydrodynamic calculations of the collapse of both singular and 
nonsingular logatropic spheres. 

In the SLS collapse solution an expansion wave moves outward into a static cloud and 
sets the gas into motion towards the central star. Both the speed of the expansion wave and 
the mass accretion rate increase with time. The radius of the expansion wave is given by 

T,^^\{27:GP,Y/H\ (1) 

where t is the time elapsed since the onset of collapse. 

Outside the radius of the expansion wave (r > Tew), the SLS envelope is static, 

^^SLS = 0, (2) 



and the density is given by 

PSLs(r) = (Po/27rG)VV-\ 



(3) 



-6- 



Inside the expansion wave (r < Tew), the infall velocity, whose direction is radial, is 
given by 

^^SLs(r) = ^(27rG'Po)'/'i^^(^), (4) 
having a free-fall behavior for r <^ Tew The density is 

2 



TT 



where x = A(2-kGPo)^^^'^ r is the similarity variable, and u(x) and a(x) arc non-dimensional 
functions. With this normalization the expansion wave is located at a; = 1. The self-similar 
variable, and the density and velocity functions arc related to those tabulated by McLaugh- 
lin & Pudritz (1997) through x = 2^^'^xmp, C({x) = 2^'^Q;Mp(a;Mp), and u{x) = 2'^^'^ump{xmp), 
where the subindex "MP" labels the McLaughlin & Pudritz (1997) solution. 

As noted by Osorio et al. (1999) the SLS collapse tends to produce massive envelopes, 
since inside the radius of the expansion wave only 3% of the mass is in the central star, while 
97% is in the collapsing envelope. 

The velocity dispersion inside the envelope due to turbulent motions is obtained assum- 
ing that Alfven waves in the cloud induce fluid motions with velocity amplitudes, 6vtuT, of the 
order of the wave speed, 6vtur — va = (dP/dpsLs)^^'^ = (-Pq/psls)^^^, where the magnitude of 
the magnetic field is B — {AnPoY^'^. random polarizations and random orientations of 
the magnetic field in the cloud, the one-dimensional velocity dispersion in the line-of-sight 
direction is given by 

1 Po^'^' 



asLs(r)= --^ ' (6) 

where the factor 1/3 comes from averaging over the solid angle the square of the magnetic 
field projections along the line of sight. 

The physical structure of the SLS collapse is characterized by the constant Pq and the 
time elapsed since the onset of collapse, t. In terms of familiar variables, this elapsed time 
is given by 

t = 7 

where M is the mass accretion rate, and M^, is the mass of the central star. Pq can be written 
as 

^° " (2ii7r3Grm2M6)V3' 

where mo = 0.0302 is the reduced (dimensionless) mass. Prom now on, we will use and 
M to characterize the dynamical model. 
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The total source of heating of the envelope, Ltot, is assumed to be the sum of the stellar 
luminosity, L^,, and the accretion luminosity, L^^c = GM^^M / R^^, where is the stellar 
radius. A value of i?* = 1 x 10^^ cm has been adopted (see Osorio et al. 1999, Hosokawa 
& Omukai 2008). The stellar luminosity is related to the stellar mass using the Schaller 
et al. (1992) evolutionary tracks. The temperature of the dust grains inside the envelope, 
T(r), is self-consistently calculated from the total luminosity using the condition of radiative 
equilibrium for outer optically thin regions of the envelope, whereas for the inner optically 
thick regions the temperature is calculated from the standard diffusion approximation (see 
details in Osorio et al. 1999). 

For temperatures > 60 K and densities > 10^ cm~^, the gas and dust are well coupled 
and are described by the same temperature (e.g., Sweitzer 1978; Doty et al. 2002). Our 
models fulfill this condition and, thus, the kinetic temperature of the gas, T^, is taken to 
coincide with the temperature obtained from the dust calculations, Tfc(r) = T(r). 

Assuming that the dust grains shield the molecules from the intense central radiation 
field, the dust destruction front, at the dust sublimation temperature (~ 1200 K), determines 
the inner radius, -Rdust, of the dust and molecular envelope. The external radius, -Rext, of the 
envelope is inferred from the observations. 



2.3. Dust Continuum Emission 

To solve the radiative transfer equation for the dust emission we have the geometry 
depicted in Figure [H Given the spherical symmetry of the envelope, the emergent specific 
intensity, 1^, for a given line of sight is only a function of the impact parameter, p, the distance 
between the line of sight and the center of the envelope. Assuming that the source function of 
the thermal dust emission is the Planck function evaluated at the dust temperature, By{T), 
we have the following set of equations that allow us to obtain the emerging intensity as a 
function of frequency and impact parameter, 

h{v,p) = h^{v) e-^ + / B,{T) e-^' dr', (9) 

Jo 

t'{u,p,s') = / k(z^)psls(s,p) ds, 

J -So 

and 

/+S0 
f^i^) PsLs{s,p) ds, 
-so 

where /bg(^) is the background intensity, «:(//) is the monochromatic dust absorption coef- 
ficient per unit mass, s is the coordinate along the line of sight (r^ = p"^ + s^), positive in 
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the direction away from the observer, and whose origin is in the plane perpendicular to the 
line of sight that contains the center of the envelope (see Fig. [1]). The integration limits 
—So and +So correspond to the edges of the envelope nearest and farthest to the observer, 
respectively, being sq = \/Rcxt ~ P^- ^^e absence of a strong background source, /bg will 
be taken as the cosmological background, /bg = -By(Tbg), where Tbg = 2.7 K. 

Following Osorio et al. (1999), the dust absorption coefficient, k(z/), at high frequencies 
(z/ > 1500 GHz) is obtained from D'Alessio (1996), who considered spherical (Mie) particles 
of graphites, silicates, iron, and water ice compounds, with the standard grain-size distribu- 
tion of the interstellar medium, N{a) oc a~^'^ (Mathis, Rumpl, & Nordsieck 1977), with a 
minimum radius of 0.005 /xm and a maximum radius of 0.3 /im. The abundance and opti- 
cal constants of the compounds were taken from Draine & Lee (1984), Draine (1987), and 
Warren (1984). The value of the absorption coefficient at 1500 GHz obtained for this dust 
composition and size distribution can be extrapolated to lower frequencies in the standard 
power-law form, resulting n{v) = 0.06(z//1500 GHz)'^ cm^ g~^ for u < 1500 GHz. The value 
of the index f3 is adopted so that it is consistent with the slope of the optically thin millimeter 
emission of the observed SED, where oc We assumed a constant dust-to-gas mass 

ratio of 0.01. Depending of the wavelength, the D'Alessio (1996) opacities can be a factor of 
up to 2-3 times lower than those of Ossenkopf & Henning (1994). Both opacity models have 
been successfully used to reproduce the SED of other high mass star forming cores (Osorio 
et al 1999, van der Tak et al 2000). 

Once the emergent intensity as a function of the impact parameter is obtained, the total 
flux density of the source at a given frequency is calculated as 

^-M = ;D^y^ U^,p)pdp, (10) 

where D is the distance of the observer to the source. Finally, the resulting model SED 
is calculated and compared with the observed SED. A grid of models is run with different 
values of the free parameters M and (or, equivalently, L^,). In this way, the model (or 
models) that are consistent with the observed SED are determined. Note that in case the 
SED is not well constrained, more than one model can be consistent with the data. Further 
details on the method are described by Osorio et al. (1999). 

2.4. Ammonia Line Emission 

In this section, we describe the calculation of the emerging intensity of the ammonia in- 
version transitions arising from a spherically symmetric envelope with the physical properties 
obtained from the dust modeling using the SLS solution. 
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The physics of the ammonia molecule is studied in detail in Townes & Slawlow (1975), 
Ho (1975), and Ho & Townes (1983). Here we only summarize some basic information 
that may be relevant to our discussion. Because ammonia is a symmetric top molecule, its 
rotational energy levels are described by two quantum numbers, namely, the total angular 
momentum, J, and its projection on the symmetry axis, K {K < J). Rotational levels with 
J = K are more populated and are called metastable, while rotational levels with J ^ K are 
short-lived and are called non-metastable. All rotational levels (except those with K = 
are split into inversion doublets. The transitions between these inversion doublets are called 
inversion transitions and occur at wavelengths of ~1.3 cm. 

As shown in Appendix \^ the excitation temperature of the ammonia inversion transi- 
tions can be obtained from a two-level model. For a given point of coordinate s along a line 
of sight of impact parameter we can estimate the excitation temperature, Tex(s,p), from 
equation (]A1|) . taking the total number density of gas molecules as n{r) = pshs{r) / {fim ^n) , 
where fjjffi — 2.3 IS the mean molecular weight in the envelope and mn is the hydrogen mass. 
In the calculation of Tex we will approximate the intensity of the local radiation field by 
the background intensity {Ir{s,p) ~ Jbg). With this approximation, equation (1A1I) gives 
a lower limit for the excitation temperature. However, as shown in Appendix [Xj for the 
physical conditions of HMCs the values obtained in this way are very close to the value of 
the kinetic temperature, which, for thermal lines, is an upper limit for the excitation tem- 
perature. Therefore, this indicates that taking = Jbg is a very good approximation for the 
calculation of the excitation temperature in HMCs. Note that, with this approximation, the 
resulting excitation temperature will become a spherically symmetric function, Tex(r). 

The spectra of the inversion transitions present a hyperfine structure due to the inter- 
action of the electric quadrupole momentum of the nitrogen nucleus with the electric field 
of the electrons. This effect splits the line into five components, the central main line, and 
two pairs of satellite lines, separated ~l-2 MHz, symmetrically placed about the main line. 
These hyperfine components are further split by weaker magnetic spin interactions. Since the 
separation between the magnetic hyperfine components is only ~10-40 kHz (0.1-0.5 km s~^), 
they are distinguishable only in very high spectral resolution observations of regions with 
very low velocity dispersion. For this reason, we neglect the magnetic hyperfine structure, 
and we restrict our analysis to only the five electric quadrupole hyperfine components. 

As shown in Appendix [B] (eq. |B13j ). the absorption coefficient of the ammonia {J,K) 
inversion transition, at a given point of coordinate s along a line of sight of impact parameter 
p (see Fig. [T]), as a function of the observed LSR velocity, V, is given by 
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gjKe 



EjK/kTrot{r) 5 



{V-V,-Vc-vs-Ls{r)s/rf/2al 



(r) 



X 



g(Trot(r)) 



(11) 



where the frequency z/q corresponds to that of the main hne of the inversion transition (given 
in Table H]), Xnh3('") is the local gas-phase ammonia abundance relative to H2, cTy(r) is the 
local dispersion of the distribution of line-of-sight velocities due to turbulent and thermal 
motions, Tj.ot{r) is the local rotational temperature, gjK is the statistical weight of the (J, K) 
rotational level (given by equation [B7] ). Ejk its rotational energy (given by equation |B8j ). 
Q(Ti.ot) is the partition function (given by equation |B12] ). Xi is the relative strength of the 
ith hyperfine component {Yl^=i = 1) ^^^1 Vi its Doppler velocity shift with respect to the 
main line (both given in Table [5]), is the LSR velocity of the ambient cloud, and the 
term fsLs(^) s/r = Vg is the line-of-sight component of the infall velocity of the envelope. In 
writing equation (ITTi) it has been assumed that n = 1.2nii^, corresponding to assume a 10% 
He abundance. Note that the excitation temperature, Tex, and the Einstein spontaneous 
emission coefficient, Aui, are those of the entire inversion transition, unsplit by hyperfine 
interactions (eq. [Al] and Table HI respectively). Due to the dipolar selection rules on the 
rotational transitions, the populations of the metastable levels are mainly determined by 
collisions; therefore. Trot, that is defined by the relative populations of the rotational levels 
is expected to be similar to the kinetic temperature. In fact. Trot is a lower limit of T^, and 
its value is different for each pair of rotational levels considered (see Appendix |B]). However, 
the agreement between Trot and T^ is expected to improve for high metastable levels, due 
to the increased excitation involved, and for high densities, where all temperatures tend to 
thermalize to T^ (e.g., Walmsley & Ungerechts 1983; Danby et al. 1988; Wilson et al. 2006). 
This is the case for G31 HMC (Osorio 2000), and hereafter we will take Trot=Tfc. 

In the above equation, the dispersion of the line-of-sight velocity distribution results from 
the contribution of the macroscopic (turbulent) and microscopic (thermal) components. 



The only free parameter in equation (JTTi) is the gas-phase ammonia abundance, Xnh3('")- 
The smallest values of the ammonia abundance reported in the literature are found in the 
coolest cores of low-mass star forming regions, with typical abundances of X^Hg = 10~^-10~^ 
(e.g., Herbst & Klemperer 1973, Ungerechts et al. 1980, Estalella et al. 1993), and some 




(12) 



where crsLs('") is given by equation and o"th(r) is given by 




(13) 



where mNHs is the mass of the ammonia molecule. 
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chemical models predict abundances as low as Xnhs = 10 ^°-10 ^ (Le Bourlot et al. 1993). 
The highest values, obtained in very hot regions of massive star formation (e.g., Millar 1997, 
Walmsley 1997, Ohishi 1997), are Xnhs = 10-^-10"^ 

Neglecting chemical effects, we consider the two simplest possibilities for the function, 
^tiHsir), that describes the gas-phase ammonia abundance inside the envelope. First, we 
consider the case in which X^ii^{r) has a constant value across the envelope. As a second 
possibility, we consider the case in which the total (solid+gas) ammonia abundance remains 
constant across the envelope but the gas-phase abundance increases towards the center be- 
cause of the temperature gradient. Variations in the gas-phase molecular abundances have 
already been explored by, e.g., van der Tak et al. (2000), and Boonman et al. (2003) to 
explain the CH3OH, HCN, and H2O emission in some regions of massive star formation. We 
assume a minimum gas-phase ammonia abundance, X^^^, in the outer, cooler parts of the 
core, where the bulk of the ammonia molecules are frozen in grain mantles. The ammonia 
molecules are released from the grain mantles to the gas phase, reaching a maximum value 
of the gas-phase ammonia abundance, X^ax, in the inner, hotter regions of the core, where 
it is assumed that all the ammonia molecules are in the gas phase. Under this simple ap- 
proximation, the gas-phase ammonia abundance as a function of the distance to the center 
of the core will be given by 



where rj{n,T) is the ratio between the solid and gas phases, obtained from the balance 
between the condensation and sublimation of molecules on dust grains following the thermal 
equilibrium equation of Sandford & Allamandola (1993) (see Appendix [U]). Two possible 
cases are considered. First, it is assumed that the ammonia molecules are directly frozen in 
grain mantles, with ri{n, T) described by equation (1040 : in this case, the ammonia molecules 
are released to the gas phase at the ammonia sublimation temperature, which is of ~60 K 
for the range of densities typical in HMCs (see Fig. [S]). As a second possibility, it is assumed 
that the ammonia molecules are trapped in water ice, being released to the gas phase only 
after sublimation of water molecules, that occurs at a temperature of ~100 K as described 
by equation fICSP (see Fig. [8]). 

As in the case of the continuum emission (§ 12. 3p . we adopt the geometry of Figure (H 
and the emergent specific intensity, 4, will be given by 




X, 



-X, 



(14) 






K(y, s,p) ds, 
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/+S0 
KiV,s,p)ds, (15) 
'SO 

where Tgx is given by equation ( lAll) . and k(V, s) is given by equation ( ITT]) . As for the dust 
calculations, in the absence of a strong background source, the background intensity will be 
taken as the cosmological background, /bg = B^(T]^g), where Tbg = 2.7 K. The geometrical 
parameters are the same as in equation (Q. 

Finally, in order to make an accurate comparison with the observations, the emergent 
intensity is convolved (in 2-D) with a Gaussian with a FWHM equal to that of the observing 
beam, and a velocity smoothing is applied to reproduce the spectral resolution of the spec- 
trometer. Note that Xmin and X^ax are the only parameters to be fitted by the molecular 
line modeling. 



3. The Data 

Numerous observations of the 031.41+0.31 region are reported in the literature. The 
main results of the continuum flux density measurements in the wavelength range ~ 10 /im 
to 1.3 cm are summarized in Tabled! Most of these observations did not have the angular 
resolution required to separate the emission of the HMC from that of the G31.41+0.31 UCHII 
region, whose emission peak is only ~ 5" apart. Thus, it is likely that in the data obtained 
with angular resolution > 5" the dust continuum emission of the HMC is contaminated by 
emission originated in the nearby UCHII. Therefore, we have considered these data as upper 
limits. In addition, we note that the flux density at 1.3 cm is likely dominated by free- free 
emission. Then, we expect the dust thermal emission at this wavelength to be well below 
the upper limit set by the observations, and therefore, in our modeling, we did not attempt 
to reproduce the emission at 1.3 cm and longer wavelengths. Finally, we considered the 7 
mm data point as a lower limit since it is likely that a fraction of the flux density is missed 
in the very high angular resolution interferometric VLA observations at this wavelength. 

Observations of several line transitions from a variety of molecular species have been 
carried out towards 031 HMC. Among them, the VLA B-array observations of the NH3 
(4,4) inversion transition (VLA project code AH0483) reported by Cesaroni et al. (1998) 
(combined with C-array data from Cesaroni et al. 1994) constitute one of the best sets of 
high angular resolution data towards 031 HMC and objects of the same kind published so 
far. These combined data have an angular resolution of 0'.'63 (corresponding to a spatial 
scale of 5000 AU at the distance of 7.9 kpc of the source), high enough to spatially resolve 
the structure of the core, revealing its extreme physical conditions and their variation inside 
the core. Thus, we used this dataset to carry out the main part of our spectral line data 
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analysis. 

Figure 9c of Cesaroni et al. (1998) shows the ammonia (4,4) spectra as a function of 
the impact parameter with respect to the center of G31 HMC. As can be seen in the figure, 
the brightness temperature is very high, reaching a value of approximately 100 K for the 
line-of-sight towards the center of the core, indicating high values of the kinetic temperature. 

The optical depth is also high, as indicated by the fact that the satellite lines, whose opacity 
is ~l/60 that of the main line, reach an intensity similar to that of the main line towards the 
center of the core. Also, lines are unusually broad, resulting in the blending of contiguous 
satellite lines, indicative of large turbulent motions. This figure also indicates that the 
physical conditions in the core change as a function of the distance to the center, as shown 
by the variation of the line intensity and ratio between main line and satellites for different 
impact parameters. 

In addition to the subarcsecond resolution ammonia data of Cesaroni et al. (1998), 
Churchwell et al. (1990) and Cesaroni et al. (1992) report a multitransitional ammonia 
study of G31 HMC carried out with the Effelsberg 100 m telescope. These observations have 
an angular resolution of 40". Although the study lacks the angular resolution to be sensitive 
to variations of the line parameters across the core, providing only the integrated ammonia 
emission of the overall core, these data are useful to further test the model predictions derived 
from the analysis of the dust and high angular resolution ammonia (4,4) line data. 

In the next section we present a simultaneous analysis of the dust continuum and am- 
monia line data that will provide a powerful diagnostic to constrain the physical conditions 
inside G31 HMC. 



4. Comparison of the Model with the Observations 

4.1. Dust Continuum Emission 

The observed SED of G31 HMC is characterized by strong millimeter emission suggest- 
ing the presence of either a hot or a dense envelope (or a combination of both). Therefore, 
"a priori" we expect the models that are able to fit the observed SED to be characterized 
by a high mass accretion rate. Bearing this in mind, we ran a grid of models covering a 
relatively wide range of values of the free parameters and M. The mass and luminosity 
of the central star have been related using the evolutionary tracks of Schaller et al. (1992) 
for a solar metallicity. We tested the same values of the mass given in the tables of Schaller 
et al. in the range 1 M© < M* < 40 M© (corresponding to 0.7 L© < < 2 x 10^ Lq), 
and for the mass accretion rate we explored the range 10~^ Mq yr~^ < M < 10~^ Mq yr~^. 
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We set the value of (3 so that the slope of the model SED was consistent with that of the 
observed SED in the millimeter wavelength range, resulting in values of (3 in the range 1-1.2. 
We set the external radius to i?ext = 30,000 AU, a value inferred from the observed images 
of the source (e.g., Cesaroni et al. 1994a, Maxia et al. 2001) and required to reproduce the 
large flux densities observed in the millimeter range. 

Because of the incompleteness of the observed SED of G31 HMC, there is more than one 
model that can fit the available observational data. In Table [2] we give the parameters of a 
set of models (corresponding to values of the stellar mass listed in the tables of Schaller et al. 
1992) that are consistent with the observed SED. The higher luminosity (and hotter envelope) 
model consistent with the observed SED has a central star with a mass of M* = 25 Mq and 
an accretion rate of M = 2.6 x 10^'^ Mq yr^^. Models with > 25 Mq produce far- and 
mid-IR flux densities that exceed the observational constraints. Lower luminosity models 
(M^K < 25 Mq) that are consistent with the SED data have colder envelopes, but they have 
a higher density (see the values of the temperature and density at r = 1000 AU listed in the 
table) in order to account for the strong millimeter emission observed in G31 HMC, resulting 
in higher values for the mass of the envelope. Note that a low value of does not imply 
that the source is a low- mass protostar since in Table [2] the models with the smaller central 
mass are also the younger sources, meaning simply that the time elapsed since the onset 
of collapse, t, is still too short to accumulate a massive central object. In fact, as can be 
seen in the table, these models are associated with the more massive envelopes and would 
eventually form the higher mass stars. 

The value of the inner radius of the envelope, -Rdust, is ~ 156 AU for model I (M* = 
25 Mq) and gets smaller as the mass of the central star decreases. We note that in all these 
models the radius of the expansion wave, rgw, is smaller than the radius of the core, -Rext! 
therefore, the outermost regions of these envelopes are static. Also, the infall accretion rates 
in these models are so high that the accretion luminosity is the most important source of 
energy. For a given value of M^, the observational uncertainties in the SED data points 
result in a variation of ~ 10% in the possible values of M that are consistent with the data. 
This results also in a formal uncertainty of the order of 10% in the derived parameters listed 
in Table [2] for a given value of M^, . 

Figure [3] shows the physical structure of the SLS envelopes corresponding to the models 
listed in Table [2l As the figure illustrates, model I has the higher values of the velocity and 
temperature, and the lower values of the density and velocity dispersion. Turbulent motions 
are the most important contributors to the local velocity dispersion at large radii of the 
envelope, reaching values of ctsls > 4 km s~^ in the outermost regions. 
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4.2. Ammonia (4,4) Line Emission 

Here we obtain the NH3 emission for the models that are consistent with the observed 
SED of G31 HMC to further test their physical properties by comparing their emission 
with the VLA observations of the NH3(4,4) inversion transition reported by Cesaroni et 
al. (1998). In order to reproduce the setup of the Cesaroni et ah observations, the 
emerging ammonia intensity as a function of the LSR velocity, V, is calculated for each 
line of sight through the core and convolved with a Gaussian beam of HPBW = 0'.'63. 
The spectra have been converted into a brightness temperature scale using the relation 
Tb = 2.1(Fj,/mJy beam~^)(HPBW/arcsec)~^. A velocity range of ~ 150 km s~^ is covered, 
and Hanning smoothing is applied to obtain an effective spectral resolution of 4.86 km s~^. 
A value of Vc = 97 A km s~^ has been adopted for the LSR velocity of the ambient cloud. 
Figure 9c of Cesaroni et al. (1998) shows the spectra obtained by averaging the observed 
ammonia emission over circular annuli around the position of the core center. Therefore, in 
order to make more accurate the comparison of the model results with the observations, the 
spectra have been calculated for the impact parameters p = 0, 5000, 10000, and 15000 AU, 
corresponding to the values sampled in this figure. 

4-2.1. Constant Gas-Phase Ammonia Abundance 

We first ran a grid of cases with different values of the gas-phase ammonia abundance, 
assumed to be constant along the envelope. We tested a wide range of values, from A^Hs = 
10^^ to Anhs = 10~^, for all the models resulting from the fit to the observed SED (discussed 
in § 14. ip . After running all these cases, we conclude that none of these models is able to 
reproduce the observed spectra. 

Figure m shows the synthetic spectra obtained for model I at different impact parameters 
and for different ammonia abundances. The spectra observed by Cesaroni et al. (1998) are 
also shown for comparison. This figure shows that model I cannot reproduce the observed 
behavior of the ammonia emission as a function of impact parameter. For a low ammonia 
abundance (left column), the main line reaches the observed brightness temperature of ~ 100 
K in the p = panel; however, the satellite lines are too weak in all the panels. Increasing 
the ammonia abundance (middle and right columns) produces an increase of the intensity 
of the satellite lines, that reach a value close to the observed one (~ 100 K) in the p = 
AU panel of the right column. Nevertheless, in this panel the main line reaches a brightness 
temperature of only ~ 40 K, which is much smaller than the observed value. This occurs 
because the main line is very optically thick (reaching an opacity of Tmain — 350) and, thus, 
the brightness temperature must be roughly equal to the kinetic temperature of the outer 
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parts of the envelope model, which are cold. Note that the brightness temperature of the 
satellite lines is higher since their opacity is 60 times lower, and they can reach the inner 
parts of the envelope where the temperatures are higher. 

Models II to XI are colder than model I (see Fig. [3]) and therefore, they are worse 
in reproducing the observed behavior of the ammonia emission. As for model I, for a low 
abundance the satellite lines are too weak in all the panels. The intensity of the satellite lines 
increases with the ammonia abundance but, as in model I, for high ammonia abundances 
the main line becomes very optically thick in the p = AU panels (rmain > 380) and only 
traces cold material from the outer envelope. Since models II to XI are colder than model I, 
the line intensities are lower and remain well below the observed values. 

Clearly, there is no hope to reproduce the observations with constant values of the 
ammonia abundance. A high optical depth is required to reproduce the intense satellite 
lines, but this high opacity prevents that the main line is formed in the inner, hottest 
regions resulting in a main line brightness temperature too low. Nevertheless, the analysis 
of constant abundance suggests the alternative of having a low abundance in the outer 
envelope and a high abundance in the inner parts. In this way, the colder, outer parts of 
the envelope would have a smaller contribution to the observed emission, that would instead 
trace inner, higher temperatures. One expects that a compromise could be reached, in which 
the ammonia abundance in the outer parts of the envelope is low enough for the main line 
not to be saturated until reaching the inner, hotter regions, where the abundance may be 
high enough for even the satellite lines to become optically thick. In addition, having a 
low abundance in the outer parts of the envelope would prevent the main line to become 
saturated in the outer panels, helping to explain the observed decrease of the intensity of 
the main line as the impact parameter increases. 



4-2.2. Variable Gas-Phase Ammonia Abundance 

The rationale for a variable ammonia abundance is provided by the calculation of the 
sublimation of molecules from dust grain mantles (see Appendix [C]) . Given the temperature 
gradient in the envelope, the outer parts have temperatures below that of sublimation of 
grain mantles and it is expected that most of the ammonia molecules are trapped there. 
The ammonia gas-phase abundance should increase drastically in the inner regions of the 
envelope, where the kinetic temperature reaches the sublimation temperature of the man- 
tles and the ammonia molecules are released to the gas phase. Therefore, inclusion of the 
process of sublimation of ammonia molecules from grain mantles should provide a more re- 
alistic description of the abundance distribution required in the calculation of the ammonia 
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emission. 

In order to try to fit the observed ammonia emission using a variable gas-phase ammo- 
nia abundance, we first ran a grid of cases with the ammonia abundance calculated using 
equation (fT^ . and with the ratio, r], between the solid and gas phases given by equation (]C4I) 
of Appendix [O (that corresponds to direct sublimation of NH3 ices). This procedure assumes 
that most of the ammonia molecules are frozen into grain mantles in the outer parts of the 
envelope, being released to the gas phase at the ammonia subhmation temperature (about 60 
K; see Appendix O and Fig. [H]). In this way, the gas-phase ammonia abundance behaves like 
a step function, with a minimum value, Xmin, in the external parts of the envelope (where the 
temperature remains below the ammonia sublimation temperature), and increasing steeply 
to its maximum value, X^^^, for points inside the radius where the ammonia sublimation 
temperature is reached. 

We have explored values for the outer gas-phase NH3 abundance in the range 10^^° < 
-^min ^ 5 X 10^^, and in the range 10~^ < Xmax ^ 10~^ for the inner gas-phase NH3 
abundance. However, none of the models resulting from the fit to the observed SED (§ 14.11) 
is able to reproduce reasonably well the NH3(4,4) observations. The main problem is that 
the maximum brightness temperatures of the main line obtained with this kind of ammonia 
abundance distribution are too low (~ 60 K), to account for the peak brightness temperatures 
of ~ 100 K observed in G31 HMC Apparently, this occurs because the main line becomes 
optically thick in regions with temperatures of the order of 60 K, the temperature at which 
the bulk of ammonia molecules that were frozen in grain mantles are released to the gas 
phase. Thus, the difficulties to fit simultaneously the main and satellite lines are similar to 
those of the case of a constant gas-phase ammonia abundance. 

As already mentioned, an alternative scenario is that ammonia molecules are mixed 
with water ice in grain mantles (as suggested, e.g., by A'Hearn et al. 1987 for comets, 
and by Brown et al. 1988 and Osorio 2000 for HMCs). In this case, ammonia molecules 
will be trapped in the grain mantles until temperatures high enough for sublimation of wa- 
ter molecules are reached. Therefore, the gas-phase ammonia abundance will be described 
by equation f[T^ . where t] corresponds now to the water sublimation and is given by equa- 
tion fICSp of Appendix O As illustrated in Figure [S], the transition from low to high gas-phase 
ammonia abundance occurs at a higher temperature (~ 100 K) than in the case of pure am- 
monia sublimation. As a result, optically thick lines will trace hotter regions of the envelope, 
resulting in higher brightness temperatures. 

As in the case of pure ammonia sublimation we ran a grid of cases with values for 
the outer gas-phase NH3 abundance in the range 10~^° < Xmin ^ 5 x 10~^, and in the 
range 10~^ < Xmax ^ 10~^ for the inner gas-phase NII3 abundance. We found that models 
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II to XI cannot fit the NH3 observations. The envelopes of the lower luminosity models 
(M^, < 15 Mq) are apparently too cold to fit the observed spectra. For these models the 
sublimation temperature of 100 K is reached at a radius < 5000 AU (see Fig. [3]), and therefore 
the spectra calculated for the outer impact parameters {p > 5000 AU) behave very much 
like a case of constant gas-phase ammonia abundance equal to Xmin (since the lines of sight 
for p > 5000 AU only cross the outer regions of the envelope, whose temperatures remain 
below the sublimation temperature). Since models II to XI could not reproduce the observed 
spectra with a constant abundance, it is not expected that these models can fit the spectra in 
the outer panels, even with a variable gas-phase ammonia abundance. In summary, models 
with < 15 Mq are too cold to reproduce the observed properties of the NH3(4,4) spectra, 
even with a variable gas-phase NH3 abundance. 

In contrast, we have been able to obtain a fairly good fit to the observed spectra using 
model I {M^, = 25 Mq). In Figure [5] (central column) we show our best fit spectra, obtained 
using this model with a value of X^i^ = 2 x 10^^ and a value of Xmax = 3 x 10^^. This set of 
parameters reproduces quite well the observed spectra. It reproduces the observed intensity 
(~ 110 K) of the main line towards the center of the core [p = 0), as well as the variation of its 
intensity as p increases. It also reproduces the observed intensity of the satellite lines in the 
p = and p = 5000 AU spectra. In the p = 10000 and p = 15000 AU spectra the observed 
intensity of the satellite lines is somewhat stronger than predicted, but its signal-to-noise 
ratio is poor and therefore we consider this discrepancy less significant. Lines are broad (note 
that the satellite lines on each side of the main line appear blended), and the model also 
reproduces quite well the observed line widths (~ 12 km s~^). Interestingly, the best fit is 
obtained with values of the abundance for the cold and hot parts of the envelope that coincide 
respectively with the typical values inferred (usually assuming a constant abundance) from 
the observations of low mass (cold) and high mass (hot) star- forming cores (see § l2.4p . giving 
additional support to the goodness of the fit. 

We also show in Figure [5] two additional cases, corresponding to a lower value of the 
abundance in the hot region of the envelope, Xmax = 3 x 10^^ (left column), and to a higher 
value of the abundance in the cold region, Xmin = 2 x 10~^ (right column), to illustrate the 
behavior of the resulting spectra when varying these parameters. As Figure illustrates, the 
satellite lines are mainly sensitive to variations in X^ax (left and middle columns), whereas 
the main line is sensitive to variations in X^in (middle and right columns). At large values 
of p {p = 10000 and 15000 AU), the spectra do not change by varying X^ax because the 
temperatures along these lines of sight are lower than 100 K and, therefore, sublimation will 
not occur. Figure [S] (left column) also shows that an abundance of X^ax = 3 x 10~^ (ten times 
lower than in the best fit) in the inner regions of the envelope is high enough to reproduce the 
observed brightness temperature of the main line, but it is too low to reproduce the satellite 
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lines. On the other hand, as the right column shows, an abundance of Xmin = 2 x 10~^ in 
the outer parts of the envelope (ten times higher than in the best fit model) is adequate to 
reproduce the brightness temperature of the satellite lines, but not for the main line, whose 
peak brightness temperature decreases to a value below the observed one. This is because 
for such a high value of Xmin the main line becomes very optically thick and the observed 
emission originates in the outer (and, therefore, colder) regions of the core. 

Since the spacing between values of listed in Table [H that correspond to those 
given by Schaller et al. (1992), is too coarse, especially in the high mass range, we explored 
additional cases obtaining the relationship between and by interpolation of the tables 
of Schaller et al. In this way, we performed a fine tuning of the mass of the central star and 
determined the range of values of the parameters that is consistent with the observational 
uncertainties of the data. The results are summarized in Table [31 In summary, the dust 
continuum SED and the spatial variation of the ammonia spectra can be explained by a 
dense envelope with a central massive star of mass M* = 20-25 Mq, undergoing an intense 
accretion at a rate M = 2-3 x 10~^ Mq yr^^, and with a total luminosity of Ltot — 2 x 10^ Lq. 
The mass of the envelope is Mgnv = 1400-1800 Mq and it has a considerable static region 
that surrounds the collapsing region. The value of the magnetic field is B = 5-6 mG. The 
outermost region has a velocity dispersion of (Tsls — 4 km s~^, dominated by turbulent 
motions. The temperature of the envelope ranges from ~ 40 K (at r = 30, 000 AU) to 
~ 1200 K (at r = 130-160 AU). The gas-phase ammonia abundance inside the envelope is 
described by a step-like function with a minimum value of 1-4 xlO^® in the outer, colder 
(< 100 K) region, and 2-4 xlO"^ in the inner, hotter (> 100 K) region. This distribution 
of the gas-phase ammonia abundance results from the sublimation of ammonia molecules 
trapped in water ice grain mantles. The minimum and maximum values of the ammonia 
abundance are similar to the typical values reported in low and high mass protostellar cores, 
respectively. 

Our results for G31 HMC can be compared with the calculations for this source made 
by other authors. Cesaroni et al. (1998) carried out an analysis of the NH3(4,4) main line 
assuming local thermodynamic equilibrium and that the line is optically thick. They claimed 
that the observed gradient in the peak brightness temperature is a good estimate of the radial 
distribution of the kinetic temperature of the gas inside the core. Since the observed peak 
brightness temperature profile of G31 HMC does not coincide with the expected kinetic 
temperature profile of a sphere, these authors concluded that neither a collapsing nor an 
expanding spherical envelope can explain the properties of the observed emission. This 
conclusion was based on a large velocity gradient (LVG) hypothesis, that implies that for 
any line of sight with p the only contribution to the peak brightness temperature (at zero 
velocity relative to the ambient cloud) arises from the plane perpendicular to the line of sight 
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that contains the center of the core (since the hne-of-sight component of the infall velocity 
of this material is strictly zero). However, if lines are broad (as is the case in G31 HMC), 
the LVG hypothesis is not valid since material from outside this plane can also contribute 
to the emission at the center of the line. Furthermore, if the envelope has a static region 
(as suggested by our analysis) the peak brightness temperature would have a significant 
contribution from the material in this static region. Therefore, this analysis does not seem 
adequate to infer the radial temperature distribution in the envelope and, consequently, the 
possibility of a collapsing spherical envelope cannot be discarded. 

Cesaroni et al. (1998) favored a geometrically thin disk seen almost face on (therefore, 
with a single temperature for each line of sight) to account for the G31 HMC NH3(4,4) 
emission. However, in order to reach a very high optical depth with such a geometrically 
thin disk a very high surface density would be required. This hypothesis was not tested 
by Cesaroni et al. (1998), and it is unclear if could correspond to a realistic scenario. 
Also, the satellite lines were not considered in the analysis, and it is also unclear how the 
observed variation of the main line to satellites ratio as a function of impact parameter can 
be reproduced. 

Therefore, we conclude that although we followed a simplified approach (a spherically 
symmetric model with a single central source of luminosity), ignoring deviations from spher- 
ical symmetry produced either by rotation (disks) or due to the intrinsic elongation of the 
cloud, our model accounts for some of the basic observed properties of G31 HMC, especially 
the highest angular resolution line observations, and constitutes a step ahead in modeling 
these structures. 



4.3. Other Ammonia Transitions 

To further test our modeling, we calculated the emerging emission for ammonia transi- 
tions other than the (4,4) using the same SLS model and ammonia abundances derived by 
fitting both the SED and the subarcsecond resolution ammonia (4,4) transition data (Ta- 
ble [3]). Unfortunately, for ammonia transitions other than the (4,4) only single-dish data 
(Churchwell et al. 1990, Cesaroni et al. 1992) are available in the literature. These observa- 
tions were obtained with the 100 m telescope, with an angular resolution of 40". Therefore, a 
convolution of the model results with a Gaussian with a FWHM=40" was performed. Since 
this angular resolution is insufficient to resolve the emission of the core, only the spectra 
towards the p = position have been calculated. A value of Vc = 97.4 km s~^ has been 
adopted for the LSR velocity of the ambient cloud. 
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In Figure [6] we show the model spectra for the (1,1), (2,2), (4,4), and (5,5) transitions 
overlaid on the spectra observed by Churchwell et al. (1990) and Cesaroni et al. (1992). 
The uncertainty in the absolute calibration of the observed spectra is estimated to be about 
30% (Churchwell et al. 1990, Cesaroni et al. 1992). As can be seen in the figure, the model 
spectrum of the (4,4) line coincides pretty well (within ~10%) with the observed one. For 
the (2,2) and (5,5) lines, the brightness temperatures of the model spectra are ~ 30% lower 
than in the observed spectra, which is still within the calibration uncertainties. The ratio of 
intensities of the main to satellite lines, which is independent of the calibration, is ~ 3.5±0.9 
for the observed (2,2) spectrum (the uncertainty in the ratio is estimated from the rms of 
the spectrum), in agreement with the value of 3.9 obtained for the model spectrum. The 
ratio is ~ 3 ± 0.8 for the observed (5,5) spectrum, which is also in agreement with the value 
of 3.5 obtained for the model spectrum. The observed line width (FWHM) of the (2,2) main 
line is ~ 6 km s~^, somewhat smaller than the value of ~ 10 km s~^ obtained for the model 
spectrum. Finally, the observed line width of the (5,5) main line is ~ 9 km s~^, which is 
very similar to the value of ~ 10 km predicted by the model. Given the coincidence of 
line ratios and line widths, and taking into account the uncertainty in the calibration, we 
consider that the model predictions are roughly consistent with the observations of the (2,2), 
(4,4), and (5,5) transitions. However, uncertainties in the calibration are large (30%) and 
a definitive conclusion is difficult to attain with the present data. For the ammonia (1,1) 
line, the brightness temperature of the model spectrum is significantly lower than that of the 
observed spectrum, and the shape of the spectrum is also different. We believe that for the 
ammonia (1,1) transition, that traces the colder gas, there is likely a significant contribution 
of cold molecular gas from outside the core (see Hatchell et al. 2000), which is not considered 
in our modeling. Also, the assumption Trot=7fc used in deriving our ammonia spectra could 
affect these results. 

We conclude that our model parameters, obtained by fitting the dust and high angular 
resolution NH3(4,4) transition data, can also reproduce, within the observational uncertain- 
ties, the intensities, line profiles, and main to satellite line ratios of the multitransitional 
ammonia data. Our modeling is carried out by fitting essentially four parameters, namely, 
the mass accretion rate and the mass of the central star (or, alternatively, the stellar luminos- 
ity) in a physically self-consistent collapse model of a SLS envelope, and the minimum and 
maximum gas-phase ammonia abundances inside this envelope. We conclude therefore that 
in massive protostars, where an important temperature gradient is expected to be present, 
an analysis in terms of a realistic, self-consistent model that takes into account the gradients 
in the physical properties inside the core is necessary to properly interpret the observational 
data. 

The modeling we carried out for G31 HMC, including the continuum SED as well as 
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the high angular resolution properties of the molecular emission, is, to our knowledge, one 
of the more complete ever made for a hot molecular core. This source appears to be a good 
laboratory to test the very early stages of massive star formation. 

5. Summary and Conclusions 

Our main results can be summarized as follows: 



1. We have succeeded in developing a self-consistent model of the physical structure of a 
hot core, described as a spherically symmetric SLS envelope infalling onto a massive 
star. The model predicts the continuum emission SED as well as the shape and intensity 
of ammonia lines. 

2. This model is able to reproduce the general observed properties of G31 HMC. In 
addition to the SED, it can reproduce the strong intensity of the main and satellite 
ammonia (4,4) lines, as well as their spatial variations as observed with subarcsecond 
angular resolution. The same model can roughly reproduce, within the observational 
uncertainties, the single-dish ammonia (2,2), (4,4), and (5,5) data. The ammonia (1,1) 
data likely have a significant contribution of cold gas outside the core that has not 
been included in the model. 

3. Although the SED fitting alone cannot constrain adequately the parameters of G31 
HMC, the simultaneous fitting of the SED and the ammonia emission allows us to 
obtain a preliminary estimate of the physical parameters of this source. Our best fit 
model for G31 HMC has a central protostar of 20-25 Mq, with an age of 3-4 xlO^ yr, a 
high mass accretion rate of 2-3 xlO~^ Mq yr~^, and a total luminosity of ~ 2 x 10^ Lq. 
The magnetic field is 5-6 mC. The mass of the envelope is large, 1400-1800 Mq, and its 
temperature is high, ranging from ~ 40 K at the outer radius of the envelope (30,0000 
AU) to ~ 1200 K at the inner radius (130-160 AU). 

A steep increase of the gas-phase ammonia abundance, from ~ 2 x 10~^ in the outer 
parts of the envelope to ~ 3 x 10~^ in the inner parts, is required to reproduce the 
observed spatial variations of the high angular resolution ammonia spectra. This is 
naturally explained as a result of the release (at a characteristic temperature of ^ 
100 K) of ammonia molecules trapped in water ice grain mantles. Therefore, taking 
into account variations in the gas-phase molecular abundances because of sublimation 
and freezing-out of molecules onto grain surfaces appears to be necessary to properly 
reproduce the highest angular resolution observations, and will probably be a key 
ingredient in the modeling of future line observations with ALMA. 
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4. Although the actual scenario for G31 HMC is certainly complex (binaries, rotation, 
outflows, ...), our simple model can provide a rough description of the main features of 
the star formation process in this source, revealing that massive protostars appear to 
be able to excite high excitation ammonia transitions, not only at circumstellar scales, 
but also in the larger scale (~ 30,000 AU) infalling envelopes. 

However, the SED of G31 HMC is still incomplete, lacking high angular resolution 
data at wavelengths shorter than 880 /im. Also, high angular resolution line data are 
only available for the ammonia (4,4) transition, and the calibration uncertainties in the 
lower angular resolution observations of other transitions are still too large to properly 
constrain the models. Higher quality data as well as an improved modeling would be 
necessary to make firm statements on the physical properties of this source. 
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A. Excitation Temperature of the Ammonia Inversion Transitions 

Since the inversion transitions inside a given rotational level of the ammonia molecule are 
more frequent than the transitions to other rotational levels, we approximate any metastable 
inversion doublet by a two-level model and we estimate the excitation temperature, Tex, of its 
inversion transition by considering a detailed balance between excitation and deexcitation, 
resulting 

(1 + Auin-^^-^) e^^^l'^^^ - 1 



ex 




(Al) 



(l/2)A„^c2/i-Vo-3n-i7-,^/,e'^-o/fcT, + i 

where z/g is the rest frequency of the inversion transition, n is the total number density of 
gas molecules (mainly H2 molecules, which dominate the collisions), is the intensity of 
the local radiation field, is the kinetic temperature, ^ui is the coUisional deexcitation rate 
coefficient, and A^i is the Einstein spontaneous emission coefficient from upper to lower level 
of the inversion doublet. 



In Table m we list the frequencies adopted for the inversion transitions of the metastable 
levels (1,1) to (6,6) of the ammonia molecule (from Pickett et al. 1998). 
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The coUisional deexcitation rate coefficient (from Ho 1977) is taken as 



^ 2.27x10-" 4^ , (A2) 



cm^ s^^ / \^ , 

for all the inversion transitions, since Sweitzer (1978) and Danby et al. (1988) show that it 
does not change significantly from one transition to another. 

The Einstein spontaneous emission coefficient from upper to lower level of the inversion 
doublet is given by 

_ 647rV3 ^ 

where /i^ is the electric dipole moment of the ammonia molecule whose value is fi^ = 1.476 
D (Cohen & Poynter 1974). The values of Aui for metastable states with J < 6 are listed in 
Table H 

As an example, Figure [7] shows the excitation temperature for the ammonia (4,4) in- 
version transition as a function of the density for a range of kinetic temperatures from 10 
to 1000 K. As can be seen from the figure, the transition is well thermalized for the density 
and temperature ranges of HMCs. A similar result is obtained for transitions (1,1) to (5,5). 



B. Calculation of the Absorption Coefficient of the NH3 Inversion Transitions 

The absorption coefficient corresponding to a transition between two levels is given by 
(Estalella & Anglada 1999): 

<V) = ^A^^n^ ^e^-^o/kn. _ 1) 0(v^), (bi) 

where is the number density of particles in the upper level, and (piV) is the profile function 
that contains the dependence on the line-of sight velocity. The profile function is actually 
the line-of-sight velocity distribution function, and is normalized so that (f)(y)dV = 1. 

For the inversion transition of the (J, K) rotational level of the ammonia molecule, the 
absorption coefficient is given by 

3 5 
f^{V) = -^A^mu (e'^'^o/^^- - 1) Y^XiUV). (B2) 

where z/q is the rest frequency of the main line {i = 3) of the inversion transition (the values 
adopted are listed in Table H] for metastable levels up to J = 6), Aui is obtained from equa- 
tion flASp and listed in Table HI riu is the number density of ammonia molecules in the upper 
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level of the (J, K) inversion doublet, and Tgx is the excitation temperature that describes 
the relative populations of the inversion doublet and is calculated using equation flAip . Due 
to the quadrupole hyperfine structure of the ammonia inversion doublet, the profile func- 
tion is given by (f){V) = Yl^=i^i4'i{V), where Xi is the LTE relative intensity of the zth 
hyperfine component given in Table [5l and (piiV) its profile function that is assumed to 
be a Gaussian. Since the profile function of the overall inversion transition is normalized, 
I-^ Yl^=i^i4'i{^)d^ = 1, and because Yl^=i^i — 1; profile function of each individual 
hyperfine component should be also normalized, so that, (j)i(y)dV = 1. Therefore, 

= ( ^y^'-e-(^-^-^-^^)V2-v, (B3) 

where Vi is the Doppler velocity corresponding to the frequency shift of each hyperfine 
component with respect to the main line (V3 = 0), Vc is the line-of-sight LSR velocity of the 
ambient cloud, Vs is the line-of-sight component of the systematic velocity field, and ay is 
the local dispersion of the distribution of line-of-sight velocities due to turbulent and thermal 
motions. 

The density of molecules in the upper level of the inversion doublet, that is required in 
equation (lB2p . can be written in terms of the density of molecules in the {J,K) rotational 
state, rijK = nu + rii, using the Boltzmann equation and the excitation temperature. 

On the other hand, the total density of ammonia molecules is given by the sum of the 
densities of all the rotational states, 

00 J' 

^NHa = X] X] njiK'- (B5) 
J'=QK'=0 

Using the Boltzmann equation, the density of ammonia molecules in the rotational state 
(J, K) can be related to the density in any other rotational state (J', K') through 

njK 9jk 

where qjk and gj'K' are the statistical weights of the two rotational states, E'jk and Ejik' 
their rotational energies, and Tjk,j'k' is the corresponding rotational temperature. The 
statistical weights are given by the relations 

4(2J+1), K^3 
gjK= { 8(2J+1), K = 3, K^O (B7) 
4(2J+1), K = 0, 
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and the rotational energies are given by 

EjK = h[BJ{J+l) + {C-B)K^] {J = 0,1,2,...-K = 0,±1,...,±J), (B8) 

where B = 2.98 x 10^^ Hz and C = 1.89 x 10^^ Hz are the rotational constants of the 
ammonia molecule (Townes & Schawlow 1975). 

Solving for rij'x' in equation flB6p and substituting in equation fIBSp one obtains 

oo J' 

^NH3 = ^ E E ^/.'K'e-(^^'-'-^-)/'=^-.^'-' . (B9) 

^''^ J'=OK'=0 



Solving the above equation for ujk and substituting the result in equation (]B4p . the 
population of molecules in the upper level of the inversion doublet n„ can be obtained in 
terms of ^nhs, 

^ nNHsdJK CBlOl 

- (1 + e^^o/kn.) Ei'=o9J'K'e-^'''''''-'''''^^'^''''''''' ■ 

This equation can be finally written in terms of the number density of H2 molecules, , 
using the fractional ammonia abundance, Xnhs = '^nh3/''^H2- Furthermore, if we assume that 
the relative populations of all the rotational states are characterized by the same rotational 
temperature, that we call Trot (this is usually the case in high density regions where collisions 
dominate and all temperatures thermalize to T^), the above equation can be written as 

where Q{T) is the partition function defined as: 

00 J' 

^(^) = E E dJ'K'C-''^'-'^'^. (B12) 

J'=0 K'=0 

In the calculation of Q{T) usually only metastable rotational states (J = K) need to be 
included, since non-met ast able states are short lived and in general are not significantly 
populated. The sum is extended up to a value of J' high enough so that the contribution of 
higher levels is negligible. 

Substituting equations ( 1B11I) and (1B3I) into ( ]B2I) . the expression for the absorption 
coefficient for the inversion transition of the (J, K) rotational state of the ammonia molecule 
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becomes 



X y ,,e-^y-y^-y^-y^?i^^i . (bi3) 



=1 



C. Condensation and Sublimation of Moleculeir Species 

Let's assume a molecular species whose total amount remains constant, while the solid 
versus gas phase ratio, rj, changes as a result of condensation and sublimation processes due 
to variations in temperature and density. Let's further assume that there is a small residual 
fraction of molecules that remain in the gas phase even at temperatures well below the 
characteristic sublimation temperature, as observations suggest that the abundance of some 
molecular species does not drop to zero in the very low temperature regions. Let's call Xmin 
the minimum gas-phase molecular abundance, reached at low temperatures (where ^ 1), 
and Xinax the maximum gas-phase abundance, reached in the hot regions = 0), where all 
the molecules are in the gas phase. Therefore, Xmoi, the gas-phase molecular abundance at 
a given point will fulfill the following equation: 

^mol + ?7(^mol ~ ^min) — ^max, (Cl) 

where the first term corresponds to the fraction of molecules in the gas-phase and the second 
term corresponds to the fraction in the solid phase, after correcting for the residual fraction of 
molecules (corresponding to Xmin) that do not follow the condensation/sublimation process. 

Therefore, the gas-phase molecular abundance can be estimated as 

Y _ ^max ~ -^min ^ (nn\ 

^mol — --, — ; r ^min- V^^) 

l + T] 

Following Sandford & AUamandola (1993), the ratio, 77, of molecules in the sohd to gas 
phase, assuming that only thermal processes are involved, can be obtained as 



V = I ) ( ZT^J eME/m) (C3) 



where Ug is the number density of dust grains, Gg is the grain radius, is the lattice vi- 
brational frequency of the molecule, mmoi is the mass of the molecule, and E is the binding 
energy of the molecule on the ice surface. For the physical conditions of HMCs the expo- 
nential term, defined by the binding energy E, dominates the behavior of this equation. 
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Assuming a given dust-to-gas mass ratio, R, the density of dust grains can be written in 
terms of the density, n, of gas molecules as Ug = SfimmuRn/ (AnagPg), where /i^ is the 
mean molecular weight, mn is the mass of the H atom, and Pg is the density of the material 
constituting the dust grains. Assuming an average radius of the grains ^ 0.1 pm (see 
§ 12. 3p . a ratio R = 1/100, and a density ~ 3 g cm~^ (Draine & Lee 1984), we obtain 
ng ~ 3 X 10~^^n and equation (]C3P can be written in terms of the gas density. 

For the case of sublimation of ammonia ices, substitution of the values of z/„ = 3.45 x 10^^ 
Hz and E/k = 3075 K (Sandford & Allamandola 1993) into equation (1C3I) gives 

= 5 X 10-3° (-^) (^^^ ^ exp(3075 K/T,). (C4) 

In the case that the ammonia molecules are trapped into water ice, we substitute the 
corresponding numerical constants for water, namely z/y = 2 x 10^^ Hz and E/k = 5070 K 
(Sandford & Allamandola 1993) into equation (1C3|) . and the equation becomes 

= 8 X 10-3° (^) (^^^ ^ exp(5070 K/T^). (C5) 
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Table 1. Compilation of the Observational Dust Data of G31 HMC 







Angular 


r lux 


Aperture 




A 




Resolution 


Density 


DIZC 




(lim) 


Instrument 


( J 


(Jy) 


( J 


Keis. 


12 


MSX 


18 


2.0 ± O.P 


~ 18 


1 


12 


IRAS 


30 


4.1 ± 0.3^ 


~ 30 


2 


21 


MSX 


18 


15.0 ± 0.9*^ 


~ 18 


1 


25 


IRAS 


30 


52 ± 4*^ 


~ 30 


2 


60 


IRAS 


60 


1090 ± 170*^ 


~ 60 


2 


100 


IRAS 


120 


2820 ± 400^ 


~ 120 


2 


450 


SCUBA 


9 


227 ± 56^ 


~ 150 


3 


850 


SCUBA 


15 


55 ± 3^ 


^ 150 


3 


880 


SMA 


0.8 


21 ± 4 




4 


1300 


OVRO 


3.6 


4.0 ±0.8 


~ 7 


5 


1300 


3m IRTF 


90 


14.2 ±0.5*^ 


~ 90 


6 


1350 


SCUBA 


22 


4.9 ± 1.0^ 


~ 22 


3 


1400 


PdBI 


0.7 


4.4 ±0.9 


~ 10 


7 


1400 


BIMA 


0.5 


3.6 ±0.7 


~ 3 


8 


2000 


SCUBA 


34 


2.9 ±0.6'' 


- 34 


3 


2700 


PdBI 


2.1 


0.7 ±0.2^ 


~ 5 


9 


3300 


PdBI 


1.8 


0.6 ±0.2'^ 


~ 10 


7 


3400 


OVRO 


4.7 


0.28 ± 0.18^= 


~ 10 


5 


7000 


VLA 


0.05 


0.0034 ± 0.0005^ 




10 


13000 


VLA 


2.2 


0.05 ± 0.01^ 


~ 5 


11 



''Taken as an upper limit, because of possible contamination from nearby 
sources due to the poor angular resolution (> 5") of the observation. 

'^Wc applied a rough correction to account for the expected free-free contribu- 
tion at millimeter wavelengths from the nearby HII region, estimated from the 
map shown in Fig. le of Cesaroni et al. 1994a. 

'^The reported value was obtained after a detailed subtraction of the expected 
free- free contribution of the nearby HII region (see Maxia et al. 2001). 

^ower limit, because a fraction of the flux density is likely resolved out (see §3). 
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*^Upper limit, because the emission at this wavelength is probably dominated 
by free-free emission from the nearby HII region. 

References. — (1) Crowther & Conti 2003; (2) IRAS PSC; (3) Hatchcll ct al. 
2000; (4) J.M. Girart, priv. comm.; (5) Maxia et al. 2001; (6) Chini et al. 1986; 
(7) Beltran et al. 2005; (8) F. Wyrowski, priv. comm.; (9) Cesaroni et al. 1994b; 
(10) S. Kurtz, priv. comm.; (11) Cesaroni et al. 1994a. 



Table 2. Parameters Derived from the SED Fitting of G31 HMC' 











/3d 










"■1000 ' 


Viooo^ 


c^iooo'^ 






' cw 


Po° 


BP 


Model 


(Mq) 


(Mq yr- 






(AU) 


(Lq) 


(Le) 


(K) 




(km s~^) 


(km s~^) 




(yr) 


(AU) 


(dyn cm~^) 


(mG) 


I 


25 


2.6 X 10" 


-3 


1.0 


156 


79000 


1.5 X 10^ 


342 


3.1 X lO'^ 


4.9 


0.78 


1.5 X 10^ 


3.9 X 10* 


2.3 X 10* 


2.2 X 10"'' 


5.2 


II 


15 


2.3 X 10" 


-3 


1.0 


115 


21000 


7.5 X 10* 


296 


4.2 X lO'^ 


3.2 


0.91 


2.1 X 10-^ 


2.7 X 10* 


1.5 X 10* 


4.0 X 10"-'' 


7.1 


III 


12 


2.0 X 10- 


-3 


1.0 


97 


10300 


5.3 X 10* 


263 


4.4 X lO'^ 


2.7 


0.95 


2.2 X 10^ 


2.4 X 10* 


1.3 X 10* 


4.6 X 10^^ 


7.6 


IV 


9 


1.8 X 10- 


-3 


1.0 


82 


4000 


3.6 X 10* 


234 


4.9 X lO'^ 


2.1 


1.03 


2.6 X 10^ 


2.0 X 10* 


1.1 X 10* 


6.1 X 10-s 


8.7 


V 


7 


1.6 X 10- 


-3 


1.0 


70 


1800 


2.5 X 10* 


212 


5.7 X lO'^ 


1.7 


1.05 


2.8 X 10^ 


1.7 X 10* 


8.8 X 10^ 


7.4 X 10-'^ 


9.7 


VI 


5 


1.5 X 10- 


-3 


1.0 


62 


550 


1.7 X 10* 


194 


7.3 X lO'^ 


1.1 


1.20 


3.6 X 10^ 


1.3 X 10* 


6.6 X 10^ 


1.22 X 10-5 


12.4 


VII 


4 


1.4 X 10- 


-3 


1.0 


54 


240 


1.2 X 10* 


180 


8.3 X lO'^ 


0.9 


1.22 


3.9 X 10^ 


1.2 X 10* 


5.7 X 10^ 


1.42 X 10-5 


13.4 


VIII 


3 


1.2 X 10- 


-3 


1.0 


46 


80 


7.9 X 10^ 


163 


9.6 X lO'^ 


0.7 


1.28 


4.4 X 10^ 


1.0 X 10* 


4.6 X 10^ 


1.81 X 10-5 


15.2 


IX 


2 


1.2 X 10" 


-3 


1.1 


43 


16 


5.2 X 10^ 


160 


1.5 X 10* 


0.3 


1.51 


6.6 X 10^ 


6.7 X 10^ 


3.1 X 10^ 


3.97 X 10-5 


22.4 


X 


1.5 


1.1 X 10- 


-3 


1.1 


37 


5 


3.5 X 10^ 


145 


1.8 X 10* 


0.2 


1.60 


7.7 X 10^ 


5.6 X 10^ 


2.5 X 10^ 


5.40 X 10-5 


26.1 


XI 


1 


1.0 X 10- 


-3 


1.2 


33 


0.7 


2.2 X 10^ 


140 


2.6 X 10* 


0.1 


1.84 


1.1 X 10* 


4.0 X 10^ 


1.7 X 10^ 


1.02 X 10-* 


35.8 



^Parameters of SLS models that are consistent with the observed SED of G31 HMC. Models were obtained adopting a given value of the mass of the central star, M, , 
and fitting the observed SED taking the mass accretion rate, M, as the only free parameter. An outer radius of the envelope Rext = 30000 AU and a radius of the central 
star Rf = 1 X 10^"^ cm have been assumed in all the models (see ij 12.21 1. For a given value of M, , the observational uncertainty in the data points results in a formal 
uncertainty of about 10% in the values of M obtained, and in the remaining derived parameters. 

''Mass of the central star. Possible values are in the range 25 Mq > M* > 1 Mq. The values listed correspond to those given in the evolutionary tracks of Schaller et 
al. (1992) in this range. 

"^Mass accretion rate. 

^Index of the power law that describes the dust absorption coefficient for u < 1500 GHz, inferred from the slope of the observed SED in the millimeter regime. 

"^Inner radius of the envelope, calculated as the radius of the dust destruction front, at the dust sublimation temperature of 1200 K (see procedure in Osorio et al. 1999). 

'uminosity of the central star obtained from the stellar mass using the evolutionary tracks of Schaller et al. (1992). 

s Accretion luminosity. Lace = GMA/»/i?». 

''Temperature at radius r = 1000 AU. 

'Number density of gas molecules at radius r = 1000 AU. 

JInfall velocity at radius r = 1000 AU. 

'^Dispersion velocity at radius r = 1000 AU. 

'Mass of the envelope, obtained by integration of the density distribution. 
""Time elapsed since the onset of collapse. 



"Radius of the expansion wave obtained from equation JTJ (see 12.21 1. 
"Pressure scale of the SLS model, obtained from equation ((Sj (see § \2. 
PMagnetic field, obtained as _B = {AttPq)'^^^ . 
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Table 3. Parameters of the Best Fit Model 
of G31 HMC^ 



£1 T m fii" f^T 


A/q 1 n 

V CtlLiC 


J?gxt 


30000 AU^ 




1 X 10^2 cjn^ 


/3 


1.0*^ 




20-25 Mq 


M 


2-3 xlO-3 Mq yr-i 


-^dust 


130-160 AU 


L* 


5-8 xlO^ Lq 


-^acc 


1.0-1.5 xlO^ Lq 




40 K 


T'(-Rdust)'^ 


1200 K 




6000-6500 AU 




1400-1800 Mq 




3-4 xlO^ yr 


''ew 


1.9-2.3 xlO^ AU 


Po 


2-3 xlO-6 dyn cm-^ 


B 


5-6 mG 


Y ■ f 


1-4x10"^ 


X § 


2-4x10-6 



^Summary of physical parameters of the SLS 
model that fits the observed SED and the high 
angular resolution NH3 (4,4) data of Cesaroni et 
al. (1998). 

'^Adopted value. 

^Temperature at the outer radius of the enve- 
lope. 

^Temperature at the inner radius of the enve- 
lope. 

^Radius where the temperature reaches a 
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value of 100 K, the sublimation temperature of 
H2O ices. 

^Ammonia abundance relative to H2 in the 
outer part of the envelope (r > rioo k)- 

^Ammonia abundance relative to H2 in the 
inner part of the envelope (r < rioo k)- 
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Table 4. Parameters of Inversion Transitions of the Ammonia Molecule 











(J,K) 


(GHz) 


(s-i) 




(1,1) 


23.6944955 


1.66838 X 10- 


-7 


(2,2) 


23.7226333 


2.23246 X 10- 


-7 


(3,3) 


23.8701292 


2.55865 X 10- 


-7 


(4,4) 


24.1394163 


2.82264 X 10- 


-7 


(5,5) 


24.5329887 


3.08642 X 10- 


-7 


(6,6) 


25.0560250 


3.38201 X 10- 


-7 



^From Pickett et al. 1998. 
^From equation flA3P 
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Table 5: Electric Quadrupolc Hypcrfine Structure of the NH3 Molecule 











(MHz) (km s-i) 


Xi Notes 


(1,1) 



1 


0^1 




1.531 


-19.37 




0.11111 


external satellite 


2 


2 -+ 1 




0.613 


-7.75 




0.13889 


internal satellite 


3 


1^1 + 


2^2 








0.08333+0.41667 


= 0.50000 


main line 


-1 


1 2 




-0.613 


7.75 




0.13889 


internal satellite 


5 


1^0 




-1.531 


19.37 




0.11111 


external satellite 


(2,2) 


1 


1 -+ 2 




2.04 


-25.78 




0.05000 


external satellite 


2 


3^2 




1.31 


-16.55 




0.05185 


internal satellite 


3 


+ 


2-»2 + 3->3 








0. 15000+0.23148+0.41481 


= 0.79629 


main line 


4 


2^3 




-1.31 


16.55 




0.05185 


internal satellite 


5 


2 -» 1 




-2.04 


25.78 




0.05000 


external satellite 


(3,3) 


1 


2^3 




2.30 


-28.88 




0.02645 


external satellite 


2 


4^3 




1.71 


-21.47 




0.02678 


internal satellite 


3 


2^2 + 


3^3 + 4^4 








0.21164+0.28009+0.40178 


= 0.89352 


main line 


4 


3^4 




-1.71 


21.47 




0.02678 


internal satellite 


5 


3^2 




-2.30 


28.88 




0.02615 


external satellite 


(1,1) 


1 


3^4 




2.45 


-30.43 




0.01620 


external satellite 


2 


5^4 




1.95 


-24.21 




0.01629 


internal satellite 


3 


3^3 + 4^4 + 5-+5 








0.24305+0.30083+0.39111 


= 0.93500 


main line 


4 


4^5 




-1.95 


24.21 




0.01629 


internal satellite 


5 


4-»3 




-2.45 


30.43 




0.01620 


external satellite 


(5,5) 


1 


4^5 




2.57 


-31.40 




0.01090 


external satellite 


2 


5^6 




2.12 


-25.91 




0.01094 


internal satellite 


3 


4^4 + 


5^5 + 6^6 








0.26181+0.31148+0.38299 


= 0.95629 


main line 


4 


6^5 




-2.12 


25.91 




0.01094 


internal satellite 


5 


5^4 




-2.57 


31.40 




0.01090 


external satellite 


(6,6) 


1 


5^6 




2.63 


-31.46 




0.00783 


external satellite 


2 


6 7 




2.25 


-26.92 




0.00785 


internal satellite 


3 


5^5 + 


6-*6 + 7->-7 








0.27422+0.31765+0.37677 


= 0.96863 


main line 


4 


7^6 




-2.25 


26.92 




0.00785 


internal satellite 


5 


6 -» 5 




-2.63 


31.46 




0.00783 


external satellite 



-40- 




Fig. 1. — Geometry of the HMC The plane containing the center of the core and the hne of 
sight is shown. i?dust ^-nd i?ext ^-re the inner and outer radii of the envelope, respectively, p 
is the impact parameter of the line of sight, r is the radius of a given point, s its coordinate 
along the line of sight, and — sq, +so are the coordinates of the envelope edges. 
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Fig. 2. — Observed flux densities of G31 HMC (see Tabled]) and predicted SED for model I 
(M, = 25 Mq; solid line) and models II-XI (M, = 15-1 M©; dashed lines). 
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Fig. 3. — Physical structure of the HMC for model I (M* = 25 Mq; solid line) and models 
II-XI {M^ = 15-1 Mq] dashed lines), (a) Dust temperature as a function of radius; (b) 
Gas density as a function of radius; (c) Infall velocity as a function of radius; (d) Turbulent 
velocity dispersion as a function of radius. 




Fig. 4. — Synthetic spectra of the NH3(4, 4) transition (sohd hne) for model I (M* = 25 Mq) 
as a function of the impact parameter, for different values of the gas-phase NH3 abundance, 
assumed constant along to the envelope. The values of the abundance are 4 x 10~^ (left 
column), 4 x 10~^ (middle column), and 4 x 10^^ (right column). The values of the impact 
parameter are indicated in the upper right corner of each panel. The observed spectra 
(adapted from Fig. 9c of Cesaroni et al. 1998) are plotted in each panel as dotted lines. The 
angular resolution is 0'.'63. 
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Fig. 5. — Synthetic spectra of the NH3(4, 4) transition (sohd hne) for model I (M^, = 25 Mq) 
as a function of the impact parameter, assuming a variable gas-phase NH3 abundance along to 
the envelope. The values of the minimum, Xmin, and maximum, Xmax, ammonia abundance 
for each case are given on the top of each column. The best fit model, with Xmin = 2 x 10~^ 
and Xmax = 3 X 10~^, is shown in the central column. The values of the impact parameter 
are indicated in the upper right corner of each panel. The observed spectra (adapted from 
Fig. 9c of Cesaroni et al. 1998) are plotted in each panel as dotted hnes. The angular 
resolution is 0''63. 
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Fig. 6.— Synthetic spectra (solid line) of the NH3(1, 1), NH3(2,2), NH3(4,4), and NH3(5,5) 
transitions towards the center of the HMC, for an angular resolution of 40". The model 
parameters and the ammonia abundances are the same as in the best model obtained by 
fitting the high angular resolution NH3(4, 4) data (central column in Fig. [5]). The spectra 
observed with the 100 m telescope (Churchwell et al. 1990, Cesaroni et al. 1992) are also 
shown (dotted lines). The calibration uncertainty in the observed spectra is 30%. 
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Fig. 7. — Excitation temperature of the ammonia (4,4) inversion transition as a function of 
the molecular gas density, for different values of the kinetic temperature, obtained from a 
two-level model. A local radiation field with the intensity of the cosmological background 
field, at 2.7 K, has been assumed. For emission lines, this gives a lower limit to the excitation 
temperature. Note that the transition is well thermahzed, Tex ^ T^, for the range of densities 
and kinetic temperatures of HMCs, which is indicated by a box. 
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Fig. 8. — Gas-phase ammonia abundance as a function of temperature, obtained from the 
balance between the condensation and subUmation of molecules on dust grains following the 
thermal equilibrium equation of Sandford and AUamandola (1993) (see text). The dotted line 
corresponds to the sublimation of pure ammonia molecules, while the solid line corresponds 
to the assumption that ammonia molecules are trapped in water ice mantles, being released 
to the gas phase after sublimation of water molecules. For this temperature range, the results 
depend weakly on the density, for which the density distribution field of G31 HMC has been 
assumed. 



